clear all; clc;close all;
%Updated: 24/05/2018

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

def0=load('gdef1.txt');
gn0=load('ggn.txt');
gnaut0=load('ggnaut.txt');
dist0 = load('mcdistA2.txt');
index_bor = find(dist0(:,1)==1);
bort0 = dist0(index_bor,1);
at0   = dist0(index_bor,2);
bt0   = dist0(index_bor,3);


na=size(a,1);
nb=size(b,1);

lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

for i=1:na;
	gnaut(i)=gnaut0(i);   
    for j=1:nb;
        def(i,j)=def0((i-1)*nb + j);   
        gn(i,j)=gn0((i-1)*nb + j);   
        gnstar(i,j)=(1-def(i,j))*gn(i,j)+(def(i,j))*gnaut(i);
    end;
end;

for i=1:na-1;
    for j=1:nb;
        gndif(i,j)=gn(i+1,j)-gn(i,j);  
        gndif(i,j)=100*gndif(i,j);
        if (def(i,j)==0);   %default region
            gndif2(i,j)=2;
        else;
            if gndif(i,j)>0.001;
                gndif2(i,j)=1;
            elseif gndif(i,j)<-0.001;
                gndif2(i,j)=-1;
            else;
                gndif2(i,j)=0;
            end;
        end;
        
    end;
end;

%%for i=1:nb;
%for j=1:na;%
%	fm(j,i)    =fm0((i-1)*na + j);        %
%	fm1(j,i)   =def(j,i)*fm(j,i)+(1-def(j,i))*(-2);        
%end;
%end;

bgrid=b;
b =-bgrid/(lambda+r);  %change of sign
gdpss=4.822;
a=exp(a);

index_a=find(at0<na);
aat0=at0(index_a);

bt0 = bt0(index_a);
bbt0=-bgrid(bt0)/(lambda+r);

aat0=a(aat0);

b=100*b/gdpss;
bbt0=100*bbt0/gdpss;

figure;
contourf(b,a(1:na-1),100*gndif2(:,1:end),3);
%contourf(b(1:end),a(1:end-1),100*gndif2(:,1:end),5);
xlabel('debt (% GDP^{ss})');AX1 = gca;box on;
ylabel('y^T');%set(AX1,'YTick',-0.8:0.05:1.2);axis(AX1, [-10 40 0.88 1.128]);
hold on;
scatter(bbt0,aat0);
